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This  paper  introduces  the  “discrete-time  realization  algorithm”  (DRA)  as  a  method  to  find  a  reduced-order, 
discrete-time  realization  of  an  infinite-order  distributed-parameter  system  such  as  a  transcendental 
impedance  function.  In  contrast  to  other  methods,  the  DRA  is  a  bounded-time  deterministic  method  that 
produces  globally  optimal  reduced-order  models.  In  the  DRA  we  use  the  sample  and  hold  framework  along 
with  the  inverse  discrete  Fourier  transform  to  closely  approximate  the  discrete-time  impulse  response. 
Next,  the  Ho-Kalman  algorithm  is  used  to  produce  a  state-space  realization  from  this  discrete-time 
impulse  response.  Two  examples  are  presented  to  demonstrate  the  DRA  using  low-order  rational- 
polynomial  transfer  functions,  where  the  DRA  solution  can  be  compared  to  known  solutions.  A  third 
example  demonstrates  the  DRA  with  a  transcendental  impedance  function  model  of  lithium  diffusion 
in  the  solid  phase  of  a  lithium-ion  battery,  showing  that  a  third-order  discrete-time  model  can  closely 
approximate  this  infinite-order  model  behavior. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Battery  cell  models  are  necessary  to  help  gain  insight  into  cell 
behavior  and  to  be  able  to  derive  battery  controls  methods  that  are 
both  effective  and  efficient.  At  the  moment,  available  models  tend 
to  fall  into  one  of  two  categories:  theoretically  developed  physics- 
based  models,  and  empirically  justified  equivalent-circuit-based 
models.  The  former  are  posed  as  coupled  sets  of  partial-differential 
equations,  which  involve  considerable  computational  resources 
to  simulate,  while  the  latter  are  generally  computable  as  sets  of 
continuous-time  ordinary  differential  or  discrete-time  ordinary 
difference  equations,  which  are  relatively  simple  to  simulate.  The 
tradeoff  is  that  the  physics  models  provide  much  more  predictive 
power  and  are  able  to  apply  even  to  unusual  operating  situations, 
while  the  empirical  models  should  not  be  used  to  extrapolate 
beyond  the  original  data  used  to  generate  their  parameters. 

The  motivation  for  research  into  reduced  order  models  is  driven 
by  a  need  for  a  model  with  predictive  properties  that  closely 
match  the  partial  differential  equation  (PDE)  model  but  with  com¬ 
putational  efficiency  similar  to  equivalent-circuit  models.  Some 
examples  of  reduced  order  models  include  Subramanian’s  single¬ 
particle  model  [1  ],  Cai’s  proper  orthogonal  decomposition  method 
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[2],  and  Smith’s  ID  transcendental  transfer  function  model  [3].  Of 
these,  we  believe  that  Smith’s  approach  has  the  greatest  promise 
for  applications  that  require  frequent  recalculation  of  the  reduced 
order  model  as  parameters  in  the  cell  change  due  to  various  aging 
mechanisms. 

Smith  derives  Laplace  (frequency)  domain  transfer  functions 
(impedance  models)  between  cell  input  current  and  solid  surface 
concentration,  overpotential,  Butler-Volmer  kinetics,  and  solid- 
electrolyte  potential  difference.  He  then  uses  two  approaches  to 
reduce  these  infinite  order  equations.  In  [4],  high-order  poles  are 
truncated  and  the  closely  placed  modes  are  combined,  taking  into 
account  the  residue  at  each  pole  (this  method  does  not  work  in  gen¬ 
eral,  but  does  for  this  example).  In  [3  ],  he  developed  a  reduced  order 
model  by  selecting  poles  and  residues  to  minimize  a  cost  function 
in  the  frequency  domain  via  a  nonlinear  optimization.  This  process 
works  well  when  it  is  done  only  once  and  can  be  supervised  by  a 
knowledgeable  engineer.  However,  we  envision  the  parameters  of 
the  physics  models  changing  over  a  cell’s  life  as  it  ages,  and  the 
need  for  unsupervised  model  order  reduction.  Nonlinear  optimiza¬ 
tion  is  not  ideal  for  this  scenario:  results  are  very  sensitive  to  initial 
conditions  chosen,  a  global  optimum  is  not  guaranteed,  the  final 
reduced-order  model  dimension  is  not  at  all  clear  except  by  trial 
and  error,  and  the  optimization  process  is  unbounded  in  time  and 
computation. 

In  this  paper,  we  introduce  a  new  method  to  automatically 
reduce  a  transfer  function  to  an  optimized  reduced-order  model. 
We  call  this  method  the  “discrete-time  realization  algorithm”  or 
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Nomenclature 

as  specific  surface  area  of  porous  electrode,  m2  m-3 

cs  concentration  of  Li  in  an  electrode  particle,  mol  m-3 

cs,avg  average  concentration  of  Li  in  an  electrode  particle, 
mol  m3 

cs,e  surface  concentration  of  Li  in  an  electrode  particle, 
mol  m3 

Cm  extended  controllability  matrix  with  m  block 

columns 

Ds  diffusion  coefficient,  m2  s-1 

F  Faraday’s  constant,  96  487  C  mol-1 

Gk  system  Markov  parameter  k 

H(s)  Laplace  domain  transfer  function 

H*(s)  Laplace  domain  transfer  function  with  a  pole  at  the 
origin  removed 

H(z)  discrete  domain  transfer  function 

H(z)  discrete  domain  residual  transfer  function 

Hd\f\  discrete  Fourier  transform 

hd[n]  approximation  to  the  continuous  time  impulse 
response  at  sample  n  (sampling  period  of  T\ ) 
hs  discrete-time  system  step  response  with  sampling 

period  of  Ts 

hstep  approximation  to  the  continuous-time  step 

response  (sampling  period  of  T\ ) 

Him  block  Flankel  matrix  with  l  block  rows  and  m  block 
columns 

jLl  reaction  current  density,  A  m-3 

Oi  extended  observability  matrix 

o\  pseudo  inverse  of  the  matrix  (9/ 

o\  upward  shift  by  one  block  row  of  the  matrix 
Ti  sampling  period  used  to  approximate  the  continu¬ 

ous  time  response,  s 
^len  duration  of  sampling,  s 

Ts  sampling  period  of  final  discrete-time  state  space 

system,  s 

res0  residue  of  the  pole  at  the  origin 

Rs  particle  radius,  m 

r  radial  dimension,  m 

t  time,  s 

U  left  singular  vectors  from  the  singular  value  decom¬ 

position 

U-i  left  singular  vectors  retained  in  the  reduced  order 
model 

U2  left  singular  vectors  discarded  from  the  reduced 
order  model 

V  right  singular  vectors  from  the  singular  value 

decomposition 

V\  right  singular  vectors  retained  in  the  reduced  order 

model 

V2  right  singular  vectors  discarded  in  the  reduced  order 
model 

E  diagonal  matrix  of  ordered  singular  values 

Ei  diagonal  matrix  of  ordered  singular  values  retained 

in  the  reduced  order  model 

E2  diagonal  matrix  of  ordered  singular  values  dis¬ 
carded  in  the  reduced  order  model 
a  singular  values  of  the  system 


DRA.  The  DRA  converts  a  nonlinear  optimization  problem  into  an 
equivalent  linear  optimization  problem,  and  robustly  solves  the 
problem  to  give  a  globally  optimal  discrete-time  reduced-order 
model.  The  method  gives  insight  into  the  “best”  reduced  order 
model  dimension  and  executes  in  deterministic  time,  suitable  for 


unsupervised  operation  in  an  embedded  battery  management  sys¬ 
tem. 

The  DRA  consists  of  approximating  the  discrete-time  impulse 
response  of  a  Laplace-domain  transfer  function  and  then  using  this 
impulse  response  with  the  Ho-Kalman  algorithm  [5]  to  produce  a 
discrete-time  state  space  realization.  We  approximate  the  discrete¬ 
time  impulse  response  assuming  a  sample  and  hold  circuit  on  the 
input.  The  FIo-Kalman  algorithm  uses  the  discrete-time  impulse 
response  to  produce  a  state  space  realization. 

There  are  different  methods  to  produce  a  discrete-time  impulse 
response  from  a  continuous  time  transfer  function.  If  the  trans¬ 
fer  function  can  be  written  as  a  rational  polynomial,  well  known 
techniques  exist  to  convert  to  a  discrete  time  system  [6].  For  tran¬ 
scendental  transfer  functions,  however,  these  techniques  do  not 
work.  It  is  possible  to  approximate  the  transcendental  transfer 
function  using  the  infinite  product  expansion  and  then  to  obtain 
a  discrete  time  impulse  response  [7,8].  For  some  transcendental 
transfer  functions,  it  is  possible  to  combine  the  impact  of  multiple 
poles  in  a  method  known  as  residue  grouping  to  derive  a  reduced 
order  system  [4].  This  approach  is  limited  to  transcendental  transfer 
functions  where  the  poles  can  be  grouped  easily.  Another  approach 
presented  in  [9]  combines  the  inverse  discrete  Fourier  transform 
with  H2 :  approximations  to  produce  the  discrete-time  impulse 
response.  In  this  paper,  we  chose  to  estimate  the  discrete-time 
impulse  response  using  a  sample  and  hold  framework  presented 
in  Section  2.1. 

The  FIo-Kalman  algorithm  gives  a  state-space  minimal  realiza¬ 
tion  from  the  discrete-time  impulse  response.  Most  applications 
use  measured  input  and  output  data  to  derive  the  impulse  response 
[10].  The  “Eigensystem  Realization  Algorithm”  (ERA)  [11]  extends 
the  Ho-Kalman  algorithm  to  deal  with  noisy  input/output  data 
[12-14].  The  Ho-Kalman  algorithm  and  ERA  have  been  used  exten¬ 
sively  for  modal  analysis  to  determine  resonant  frequencies  in 
flexible  structures  in  cases  where  it  is  difficult  to  apply  an  impulse 
input  to  the  system  [  1 5-20].  Unlike  these  cases,  we  are  starting  with 
a  known  transcendental  transfer  function  and  then  estimating  the 
discrete-time  impulse  response.  The  Ho-Kalman  algorithm  is  then 
used,  not  with  measured  system  data,  but  with  the  approxima¬ 
tion  of  the  discrete-time  impulse  response  to  find  a  discrete-time 
state-space  reduced-order  model  of  the  system  described  by  the 
transcendental  transfer  function. 

We  note  that  the  DRA  has  very  general  application  to  any 
problem  that  requires  a  reduced-order  approximate  model  to  a 
higher-order  transfer  function,  and  is  ideal  for  reducing  the  order 
of  distributed  parameter  (infinite  order)  systems.  Modeling  bat¬ 
tery  dynamics  is  only  one  possible  application.  In  this  paper, 
our  concern  is  to  introduce  and  derive  the  DRA,  and  ultimately 
to  show  how  it  can  be  applied  to  one  equation  that  is  part  of 
a  battery  model.  A  future  paper  will  show  how  a  full  coupled 
reduced-order  battery  model  can  be  created  with  the  aid  of  the 
DRA. 

The  paper  is  organized  as  follows.  In  the  next  section,  we 
introduce  the  DRA  and  the  procedural  steps  to  derive  the  discrete¬ 
time  realization  from  a  Laplace-domain  transfer  function.  The 
algorithm  is  illustrated  on  a  simple  second-order  system  with  a 
continuous-time  rational-polynomial  transfer  function  and  then 
on  a  third-order  system  with  an  integrator  term.  In  both  cases, 
the  results  are  compared  to  the  exact  output  of  the  continuous 
to  discrete  conversion  using  the  zero-order  hold  method.  The 
third  example  illustrates  the  performance  of  the  algorithm  on 
a  transcendental  transfer  function  that  models  the  diffusion  of 
lithium  in  a  porous  electrode  [4].  We  demonstrate  two  methods 
to  deal  with  the  pole  at  the  origin  of  this  transfer  function.  The 
discrete-time  realization  from  the  DRA  is  compared  to  numeric 
solution  using  a  1-D  parabolic-elliptic  partial  differential  equation 
solver. 
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2.  Derivation  of  the  discrete-time  realization  algorithm 


2.  I .  Sample  and  hold  framework 


Given  a  continuous-time  transfer  function  in  the  Laplace 
domain,  H(s)  =  Y(s)/l/(s),  and  a  sampling  period,  Ts,  we  want  to 
derive  a  reduced-order  discrete-time  state-space  realization  of  the 
form 

x[k+  1]  =Ax[k]  +  Bu[k]  m 

y[k]  =  Cx[k]  +  Du[k],  1  } 

where  the  first  equation  is  known  as  the  “state  equation”  and 
describes  the  dynamics  of  the  system,  and  the  second  equation 
is  known  as  the  “output  equation”  and  describes  how  the  output 
y[k]  g  JZq  is  computed  as  linear  combination  of  the  states  x[k]  e  Kp 
and  the  input  u[k)  g  JZm .  The  matrices  A  g  npxp,B  e  npxm,C  g  JZqxp , 
and  D  g  jzqxm  are  constant  in  this  work.1  A  sufficient  condition  for 
the  DRA  to  operate  is  that  H(s)  be  an  element  of  the  Hardy  space 
Hoo,  which  implies  that  it  is  a  strictly  stable  and  proper  system.  This 
is  not  a  necessary  condition,  however,  as  we  will  later  generalize 
the  method  to  work  with  systems  having  isolated  pole(s)  on  the 
imaginary  axis.  Note  that  we  do  not  restrict  H(s)  to  be  formulated 
as  a  quotient  of  polynomials  in  the  Laplace  variable  “s”  (for  which 
well  known  methods  exist  to  find  the  discrete-time  system). 

At  the  heart  of  the  DRA  is  the  Ho-Kalman  algorithm  [5],  which 
takes  a  system’s  set  of  Markov  parameters  (which,  for  a  single-input 
single-output  system  are  equivalent  to  its  discrete-time  impulse 
response  sequence)  and  computes  from  them  a  discrete-time  state- 
space  model  of  the  system.  To  use  this  algorithm,  we  must  then  first 
compute  the  discrete-time  impulse  response  from  the  continuous¬ 
time  transfer  function.  The  discrete-time  impulse  response  itself 
is  computed  from  the  continuous-time  impulse  response,  which 
is  approximated  as  an  inverse  frequency  transform  of  the  initial 
transfer  function.  We  describe  the  algorithm  in  four  steps,  which 
we  preview  here,  and  discuss  in  more  detail  in  the  following  sub¬ 
sections. 


1  Sample  the  continuous-time  transfer  function  in  the  frequency 
domain  at  a  high  rate,  and  take  the  inverse  discrete  Fourier  trans¬ 
form  (IDFT)  to  get  an  approximation  to  the  continuous  time 
impulse  response. 

2  Compute  the  discrete-time  impulse  response  values  from  the 
continuous-time  response,  assuming  a  sample  and  hold  circuit 
connected  to  the  system  input. 

3  Generate  a  discrete-time  state-space  realization  using  the  deter¬ 
ministic  Ho-Kalman  algorithm.  This  algorithm  returns  the 
reduced  order  A,  B,  and  C  matrices  from  the  discrete-time  impulse 
response  sequence  in  Step  2.  The  order  of  the  system  is  deter¬ 
mined  from  the  ordered  singular  values  of  the  Hankel  matrix 
computed  as  part  of  the  algorithm.  The  D  matrix  is  found  by  the 
initial  value  theorem. 

4  Transform  the  state  space  system  into  the  desired  final  form  using 
a  similarity  transformation,  if  required. 


Knowing  that  we  require  a  discrete-time  impulse  response  in 
order  to  invoke  the  Ho-Kalman  method,  we  take  advantage  of  a 
well  known  result  from  discrete-time  system  theory.  Namely,  the 
discrete-time  transfer  function  H(z)  corresponding  to  a  continu¬ 
ous  time  transfer  function,  H(s),  assuming  that  the  input  to  H(s)  is 
piecewise  constant  with  period  Ts,  is  [6]: 

H(z)  =Z 

=  Z 

=  z 

This  equation  is  typically  used  for  analytic  computations,  but  we 
will  use  it  numerically.  Note  that  the  very  compact  notation  Z[  •  ] 
means:  “find  the  continuous-time  time-domain  function  corre¬ 
sponding  to  the  Laplace-transform  frequency-domain  argument, 
then  sample  the  continuous-time  time-domain  function,  then  take 
the  z-transform  of  the  samples.”  Noting  that  H(s)/s  is  the  step 
response  of  H(s),  we  see  that  Eq.  (2)  shows  that  the  discrete-time 
impulse  response  corresponding  to  H(z)  is  a  one-step  difference 
of  the  z-transform  of  the  sampled  continuous-time  step  response. 
This  implies  that 


1  -  e~sTs 


H{s) 


(1  _  e-s:rs)^ 


H(s) 


-Z~^z 


H(S) 


(2) 


h[k]  =  hs[k]-hs[k- 1],  (3) 

where  hs[k]  is  the  discrete-time  step  response  of  the  system  with 
Hs(z)  =  Z[H(s)/s]. 

We  are  interested  in  computing  h[k]  of  Eq.  (3),  so  we  must  first 
compute  hs  [/<].  This  could  be  done  by  working  with  H(s)/s,  but  doing 
so  has  numeric  issues  as  s  ->  0.  Instead,  we  recognize  that  H(s)/s 
represents  the  step  response  of  the  system,  which  can  also  be  com¬ 
puted  in  the  time  domain  by  integrating  the  impulse  response.  That 
is  the  approach  we  take  here. 

We  approximate  the  continuous-time  impulse  response  via  a 
“discrete  equivalent”  or  frequency-domain  emulation  approach 
[21  ].  This  allows  us  to  write 

Hd(z)«H(s)|  2Z_, 

's=TfIT T 


where  T\  is  an  emulation  sampling  period  (different  from  and  gen¬ 
erally  significantly  shorter  than  the  final  system  sampling  period 

Ts)-2 

We  now  recognize  that  the  discrete  Fourier  transform  (DFT)  of 
a  sequence  is  related  to  its  z-transform  via  the  relationship  [23] 

wdL n  = 


H  21 


expO‘2jr//JV) 
giZirflN  _  1 

ej2jr//N  +  ] 


0  <f<N, 


(4) 


We  note  that  a  system  having  a  pole  at  the  origin  does  not  meet  the 
strictly  stable  requirement.  However,  we  also  show  that  this  pole 
can  be  automatically  accounted  for. 

We  now  proceed  to  discuss  the  sample-and-hold  framework, 
the  Ho-Kalman  method,  and  accounting  for  a  possible  pole  at  s  =  0 
in  more  detail. 


1  Note  that  the  true  system  being  approximated  has  state  dimension  n,  but  we  are 
deriving  a  reduced  order  approximation  to  the  true  system  having  state  dimension 
p<n. 


where  N  is  the  number  of  points  chosen  for  the  underlying 
sequence,  and  is  usually  chosen  to  be  a  power  of  2  for  efficient 
computations.  Criteria  for  choosing  N  are  given  in  Section  2.2.  The 
inverse  DFT  of  Hd\f\  gives  hd[n],  which  is  the  approximation  of 


2  In  order  to  arrive  at  an  accurate  estimation  of  the  continuous  time  transfer  func¬ 
tion,  the  sampling  frequency,  F\  =  1/Ti,  must  be  high  enough  to  capture  the  system 
dynamics.  As  a  rule  of  thumb,  the  sampling  frequency  must  be  at  least  20  times  the 
as  great  as  the  bandwidth  of  the  system  to  get  an  rough  approximation  in  the  fre¬ 
quency  domain  [22].  A  higher  emulation  sampling  frequency  gives  more  accurate 
results. 
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the  continuous-time  impulse  response  at  the  emulation  sampling 
period,  T\ 


N-[ 

hd[n]  =  ^HdLne/w'. 


(5) 


/= 0 


The  cumulative  summation  of  this  impulse  response  yields  the 
approximation  to  the  continuous-time  step  response 


k- 1 


frstep[fc]  —  T]  'y  JldU]' 


(6) 


i=0 


This  result  is  interpolated  with  sample  period  Ts  to  give  hs[k].  The 
difference  between  hs[k]  and  hs[k-  1]  gives  the  impulse  response 
of  the  system  via  Eq.  (3). 


2.2.  The  Ho-Kalman  method 

The  following  discussion  of  the  Ho-Kalman  algorithm  closely 
follows  Ref.  [24].  Given  the  impulse  response  of  the  system  we  find 
the  A,  B,  and  C  state  space  matrices  of  Eq.  (1)  and  an  approximation 
p  to  the  system  order  n. 

The  Markov  parameters,  Gk ,  of  the  system  are  given  by, 


Gk  = 


r  d, 

\  CAk~h 


k  =  0 
k  =  1,2, . 


For  a  single-input  single-output  system,  these  are  simply  equal  to 
the  discrete-time  impulse-response  sequence  of  the  system  (i.e., 
Gk  =  h[k]  from  Eq.  (3)).  Multi-input  multi-output  systems  can  be 
accommodated  in  a  straightforward  way  by  both  the  Ho-Kalman 
method  and  therefore  the  DRA,  but  we  omit  the  extra  complexity 
from  our  discussion  by  focusing  on  the  simpler  case. 

Given  the  Markov  parameters,  we  can  form  the  system  block 
Hankel  matrix 


'Ci 

C2 

C3  . 

•  •  Gm 

c2 

C3 

C4  • 

•  •  Gm+ 1 

«l,m  = 

C3 

c4 

c5  • 

•  •  Gm+ 2 

C, 

C/+i 

C/+ 2 

Q+m-1 

We  see  now  that  the  required  length  of  the  discrete-time  impulse 
response  is  N=l  +  m.  Further,  a  requirement  of  the  Ho-Kalman 
method  is  that  both  m  and  n  be  at  least  as  large  as  the  reduced- 
order-model’s  order  p.  We  usually  choose  l  =  m,  so  N>2n.3  When 


where  the  extended  observability  matrix,  Oh  and  the  extended  con¬ 
trollability  matrices,  Cm,  are  defined  as, 


Oi  = 


C 

CA 

CA2 


CA‘ 


I- 1 


Cm=[B  AB  A2B  ■■■ 

If  we  are  able  to  factor  Hkm  into  Oi  and  Cm,  then  we  can  directly 
extract  the  desired  C  matrix  as  the  top  block  row  of  (9/  and  the 
desired  B  matrix  as  the  left  block  column  of  Cm.  We  can  do  further 
processing  to  find  the  A  matrix. 

To  do  the  factoring,  we  rely  on  the  singular  value  decomposi¬ 
tion  (SVD)  [25],  which  allows  us  to  write  Hkm  =  H£Vr,  where  U 
and  V  are  orthogonal  matrices,  and  £  is  a  diagonal  matrix  with 
non-negative  entries.  The  diagonal  of  £  holds  the  so-called  singu¬ 
lar  values  of  Hkm,  which  are  ordered  such  that  o\  >  a2  > . . .  >  am. 
The  magnitude  of  the  singular  values  is  a  measure  of  the  relative 
importance  of  its  corresponding  state  to  the  overall  system  behav¬ 
ior.  By  counting  the  “large”  singular  values,  we  arrive  at  an  estimate 
p  of  the  system  order.  We  then  rethink  the  singular  value  decom¬ 
position  as 


Ei  0 
0  £2 
=  i/1e1v[  +  &2e2v’J, 


W/,m  =  U,  U2 


4 


(7) 


retained 


discarded 


where  £i  e  Upxp  and  we  retain  only  Hkm  «  Hi  £i  Vf.  The  approx¬ 
imation  is  exact  if  £2  =  0*  We  know  from  the  Schmidt-Mirsky 
theorem  that  this  approach  yields  a  globally  optimum  p-rank 
approximation  [26]. 

It  is  possible  to  show  that  the  original  system’s  extended  observ¬ 
ability  matrix  and  the  extended  controllability  matrix  are  related 
to  the  singular  value  decomposition  via 


O,  =  US1/2!,  and  Cm  =  T-1  ?}I2VT , 


(8) 


where  T is  an  invertible  similarity  transformation  matrix  that  deter¬ 
mines  the  basis  in  which  the  A,  B,  C,  and  D  matrices  are  defined.  This 
basis  is  not  important  for  the  function  of  the  reduced-order  model, 
so  we  simply  set  T=/and  compute  (9/  =  Hi  £]^2  and  Cm  =  E-^V7. 

At  this  point  the  order  of  the  reduced-order  model,  p,  has  been 
determined  from  the  singular  values.  The  B  matrix  has  been  found 
by  taking  the  first  block  row  of  Cm.  Likewise,  the  C  matrix  has  been 
found  by  taking  the  first  block  column  of  (9/.  It  remains  to  find  A 
and  D. 

The  A  matrix  can  be  found  via  the  relationship, 


executing  the  Ho-Kalman  method  the  first  time,  when  the  actual 

"c 

CA 

system  order  is  uncertain,  an  assumed  upper  bound  of  the  sys¬ 

CA 

CA2 

tem  order  is  chosen  a  priori.  This  bound  can  later  be  verified  by 

CA 2 

A  = 

examining  the  singular  values  of  the  Hankel  matrix,  as  we  will 

demonstrate. 

_  CAl~2  _ 

_CAl~\ 

An  interesting  feature  of  the  block  Hankel  matrix,  which  is  crit- 

ical  for  this  development,  is  that  it  can  be  written  as  Hlm  =  Ofim.  The  matrjx  Qn  thg  left  side  of  the  equation  is  0(1  and  the  matrix  on 


3  For  transcendental  transfer  functions  having  an  infinite  number  of  poles,  we 

cannot  choose  N  to  be  twice  the  number  of  actual  poles.  We  find  that  best  accuracy 

is  achieved  when  N  is  chosen  such  that  the  approximate  discrete-time  impulse  used 

as  input  to  the  Ho-Kalman  algorithm  response  settles  to  near  its  steady-state  value. 


the  right  hand  side  of  the  equation  is  o\  where  the  f  represents  an 
upward  shift  of  the  matrix.  The  A  matrix  may  be  solved  from  this 
relationship  by  matrix  manipulations  and  back  substitution,  and 
the  solution  may  be  written  as 


71  =  0j  .a 


nt 

l-\Ul  ’ 
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where  the  symbol  f  represents  the  matrix  pseudo  inverse  [25]. 
Finally,  the  D  matrix  is  the  first  Markov  parameter  but  can  often 
be  better  found  numerically  by  using  the  initial  value  theorem  for 
a  discrete  time  system, 


2.4.  Summary  of  the  DRA 

Here,  we  put  together  all  the  steps  of  the  DRA  in  one  place, 
summarizing  the  method. 


D  =  G0  =  h[  0]  =  limH(z).  (9) 

Z — >oo 

The  singular  values  of  the  system  Hankel  order  provide  clear  insight 
into  the  relative  importance  of  each  state.  Instead  of  merely  trun¬ 
cating  all  higher  order  states  balanced  residualization  [27]  could  be 
used  to  include  some  of  the  behavior  of  these  less  important  states 
in  the  final  realization.  For  example  if  a  final  third  order  system  is 
desired,  the  Ho-Kalman  method  could  be  used  to  produce  a  sixth- 
order  system.  The  balanced  residualization  technique  could  then  be 
used  to  reduced  the  system  to  third  order,  incorporating  behavior 
from  the  least  important  three  states. 

2.3.  Dealing  with  one  or  more  poles  in  H(s)  at  the  origin 

If  the  original  system,  H(s),  has  a  pole  at  the  origin,  it  is  not 
strictly  stable,  so  violates  the  necessary  conditions  that  make  the 
DRA  work.  However,  it  is  quite  simple  to  deal  with  this  case.  We 
first  subtract  the  pole  at  the  origin  from  the  transfer  function, 
then  execute  the  DRA  on  the  residual  system,  then  compute  a  final 
discrete-time  state-space  model  that  augments  the  DRA  result  with 
additional  dynamics  to  implement  the  function  of  the  s-domain 
pole  at  the  origin  (we  apply  the  same  general  approach  to  remove 
multiple  poles  at  the  origin,  or  poles  elsewhere  on  the  imaginary 
axis  of  the  s-plane). 

A  pole  at  the  origin  is  removed  by  first  calculating  the  residue  of 
this  pole  and  then  subtracting  it  from  the  original  transfer  function. 

H*(s)  =  H(s)-^  (10) 


•  If  the  original  system  H(s)  has  a  pole  at  the  origin,  determine  res0 
using  Eq.  (11)  and  compute  H*(s)  using  Eq.  (10).  Computer  code 
can  automatically  determine  whether  H(s)  has  a  pole  at  the  origin 
by  substituting  small  values  of  s.  If  |H(s)|^  oo,  then  this  “zeroth” 
step  must  be  performed.  Perform  the  following  four  steps  of  the 
DRA  using  H*(s)  instead  of  H(s). 

•  Step  1:  Select  an  emulation  sampling  period  T i  such  that  1/Ti 
is  significantly  greater  than  the  bandwidth  of  H(s),  and  N  is  at 
least  twice  as  large  as  the  assumed  reduced-order  model  dimen¬ 
sion.  Compute  Hd[f]  using  Eq.  (4),  and  then  the  approximate 
continuous-time  impulse  response  via  Eq.  (5). 

•  Step  2:  Find  the  approximate  continuous-time  step  response  of 
H(s)  via  Eq.  (6).  Then,  interpolate  at  the  final  desired  sample  rate 
Ts  to  give  hs[k].  Difference  hs[k]  according  to  Eq.  (3)  to  yield  h[k\, 
the  discrete-time  impulse  response  for  H(s). 

•  Step  3:  Populate  the  Hankel  matrix  H^m  with  the  discrete-time 
impulse  response  values  h[k\.  Perform  the  singular  value  decom¬ 
position  of  Him  and  determine  from  the  singular  values  the 
system  order  p.  Partition  the  singular  value  decomposition  to  give 
the  matrices  Ui ,  Ei,  and  V\.  From  these,  compute  Oi  and  Cm  via 
Eq.  (8).  Compute  the  reduced-order-model  matrices  A,  B ,  and  C 
from  Oi  and  Cm.  Compute  D  via  Eq.  (9). 

•  Step  4:  Depending  on  the  application  it  may  be  beneficial  to  do  a 
similarity  transformation  of  the  state  space  system.  For  example, 
diagonalizing  the  A  matrix  could  speed  up  the  calculations,  which 
could  be  important  if  this  is  used  in  a  real  time  control  application 
[28]. 

•  Finally,  if  the  original  system  H(s)  has  a  pole  at  the  origin,  find  the 
final  state-space  form  via  Eqs.  (12)  and  (13). 


where  reso  is  calculated  as 

reso  =  lim  sH(s).  (11) 

s^-0 

The  remainder  of  the  DRA  is  executed  using  H*(s)  instead  of  H(s). 

To  see  how  to  re-incorporate  the  effect  of  the  s-plane  pole  at 
the  origin  into  the  final  reduced-order  model,  recall  that  a  pole  at 
the  origin  of  the  s-plane  corresponds  to  an  integrator.  The  discrete¬ 
time  equivalent  of  an  integrator  with  gain  res0  can  be  implemented 
as 


Xj[k  + 1]  =Xi[k]  +  Tsu[k]. 

Yi[k]  =  res0x,[k] 

Therefore,  we  can  implement  a  reduced-order  approximation  to 
the  original  system  via  the  augmented  discrete-time  state-space 
model: 


Xi[k+  1] 

1  :  0 

x,[k] 

1 

"  T,  ' 

x[k  +  1] 

0  j  A 

x[k\ 

T 

B 

*aug  [k+U  Aiug  ^aug  #aug 


u[k] 


y[k]  = 


res0 


C 


Xi\k\ 


x[k] 


+Du[k]. 


•^aug  M 


(12) 


(13) 


The  A,  B ,  and  C  matrices  in  the  above  equations  are  those  generated 
by  Step  3  of  the  DRA. 


Note  that  the  steps  of  the  DRA  require  only  standard  linear-algebra 
and  signal-processing  routines.  The  method  finds  the  globally  opti¬ 
mal  reduced-order  discrete-time  approximation  to  the  original 
continuous-time  system.  It  does  not  require  nonlinear  optimiza¬ 
tion,  and  does  not  require  iteration. 

3.  Examples  to  illustrate  the  method 

We  now  look  at  three  different  examples  to  illustrate  the  oper¬ 
ation  of  the  DRA.  The  first  two  examples  are  rational-polynomial 
transfer  functions,  which  we  use  because  we  can  calculate  the  exact 
solution  using  other  methods  [29].  We  can  then  compare  the  exact 
solutions  to  the  approximate  solutions  obtained  by  the  DRA.  The 
third  example  does  not  have  a  closed-form  exact  solution,  but  we 
can  use  a  1-D  parabolic-elliptic  partial  differential  equation  solver 
to  find  an  accurate  near-exact  solution  against  which  to  compare 
the  DRA  solution.  We  find  excellent  agreement  between  the  exact 
solutions  and  DRA  solutions  in  all  cases. 


3.1.  Example  1:  rational  polynomial  transfer  function 


The  DRA  method  is  first  applied  to  a  simple  second-order  sys¬ 
tem.  We  require  a  discrete-time  realization  with  the  a  sampling 
period  of  Ts  =  0.1  s  from  the  continuous-time  transfer  function 


Hi  00  = 


s2  +  20s  +  100 
s2  +  2s  +  8 


(14) 


This  system  has  complex  poles  at  -1  ±j2.65  rad  s-1  and  two  zeros 
at  1 0  rad  s-1 .  The  magnitude  response  of  H i  (s)  is  shown  in  Fig.  1 . 

Stepl :  The  sampling  frequency  is  selected  as  256  Hz  which  is  sig¬ 
nificantly  greater  than  the  system  bandwidth.  The  sampling  length 
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Bode  magnitude  plot  for  H  (s)  Approximate  discrete-time  impulse  resp.  for  H^s) 


Fig.  1.  Bode  magnitude  plot  for  Example  1. 


Continuous-time  impulse  responses  for  H^s) 


N  is  set  to  64  which  allows  up  to  a  32  x  32  system  Hankel  matrix 
in  Step  3.  The  transfer  function  is  sampled  at  discrete  frequencies 
according  to  Eq.  (4).  The  inverse  DFT  yields  an  approximation  to 
the  continuous  time  impulse  response.  Fig.  2  compares  the  approx¬ 
imate  continuous-time  impulse  response  computed  via  the  inverse 
DFT  to  the  exact  continuous-time  impulse  response  of  Eq.  (14).  We 
see  that  the  two  solutions  are  coincident. 

Step  2 :  The  approximation  to  the  continuous-time  step  response 
is  found  by  doing  a  cumulative  summation  of  the  impulse  response. 
The  results  for  this  example  are  shown  in  Fig.  3  and  show  excellent 
agreement  with  the  exact  step  response  of  the  continuous  time 


system.  We  linearly  interpolate  the  step  response  to  give  values  at 
integer  multiples  of  0.1  s.  The  down-sampled  step  response  is  dif¬ 
ferenced  to  yield  the  discrete-time  impulse  response  of  Fig.  4.  Again, 
there  is  excellent  agreement  between  the  approximate  impulse 
response  and  the  exact  solution,  with  the  exception  of  the  single 
point  at  t  =  0.  This  is  often  the  case  because  of  some  properties  of  the 
inverse  DFT,  but  it  causes  no  problems  since  the  impulse  response 
value  at  t= 0  is  not  used  by  the  Ho-Kalman  algorithm  in  Step  3,  and 
the  D  matrix  (which  is  equal  to  the  impulse  response  value  at  t= 0) 
is  computed  differently,  using  Eq.  (9). 

Step  3.  The  deterministic  Ho-Kalman  algorithm  is  used  to 
find  a  state-space  realization  from  the  approximate  discrete¬ 
time  impulse  response  from  Step  2.  The  approximation  of  the 
continuous-time  impulse  response  was  truncated  in  Step  1  to  64 
points,  which  allows  a  maximum  Hankel  matrix  of 32  x  32.  The  SVD 
of  the  Hankel  matrix  gives  insight  into  the  order  of  the  system.  A  log 
plot  of  the  singular  values  is  shown  in  Fig.  5.  The  first  two  singular 
values  are  almost  three  orders  of  magnitude  greater  than  the  third 
singular  value,  so  we  select  a  reduced-order  model  dimension  p  =  2. 

Truncating  to  the  first  two  states  only,  the  Ho-Kalman  algo¬ 
rithm  gives  a  state-space  realization  with  the  following  A,  B,  and  C 
matrices 

„  [  0.8656  -0.2367" 

0.8811 


The  D  matrix  is  found  from  the  initial  value  theorem  and,  for  this 
example,  is  D  =  [1]. 


a  = 


B  = 


r  — 


0.2367 

-1.624' 

0.7694 

1 


Continuous-time  step  responses  for  H^s) 


Singular  values  of  Hankel  matrix  for  H^s) 
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Fig.  3.  Step  responses  for  Example  1. 


Fig.  5.  Singular  values  indicate  only  two  significant  states  in  Example  1. 
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Final  model  discrete-time  impulse  response  for  H^s) 


Fig.  6.  Comparison  of  final  results  from  the  DRA  to  the  true  discrete-time  solution 
for  Example  1. 


The  residual  error  between  the  exact  solution  and  the  DRA  is 
seen  by  comparing  the  pole-residue  representation  of  the  discrete 
time  systems.  The  transfer  function  of  the  true  discrete-time  system 
can  be  written  as, 


Hi(z)  = 


a 


z-P\  z~P2 


where  p\  and  p2  are  poles  and  a  and  b  are  the  corresponding 
residues.  In  this  example  the  poles,  p\  and  p2 ,  can  be  solved  for  ana¬ 
lytically,  and  are  0.87335  ±j0.2366  and  the  residues  likewise  can 
be  found  to  be  1.019  =j=jl. 205.  For  the  DRA  realization,  the  poles 
are  computed  to  be  at  the  same  location  (up  to  the  fifth  signifi¬ 
cant  figure)  but  the  residues  are  computed  to  be  1.023  =pjl. 197. 
The  residuals  transfer  function  is, 


Step  4:  We  have  chosen  not  to  transform  the  system  represen¬ 
tation  in  this  example. 

A  comparison  of  the  impulse  response  of  the  DRA  discrete-time 
state-space  realization  to  the  exact  impulse  response  is  shown  in 
Fig.  6.  The  results  agree  very  well  (note  that  the  impulse  response 
value  at  t  =  0  has  been  corrected  by  the  correct  calculation  of  the  D 
matrix  in  Step  3).  Because  the  impulse  responses  agree  very  well, 
the  response  of  the  reduced-order  model  will  also  agree  well  with 
the  exact  response  for  any  input  signal  u[k]. 


Hi  (z)  = 


a  adra  _j_  ^  ^dra 
Z-Pi  Z  —  p2 


where  adra  and  bdra  are  the  residues  of  the  DRA  realization.  At 
the  256 FIz  sampling  frequency,  the  residue  errors  for  H^(z)  are 
-0.00418  =fj0.00762.  The  residual  transfer  function  is  determin¬ 
istic,  but  small  because  the  residues  of  the  true  system  and  the 
DRA  are  close.  The  error  can  be  decreased  further  by  increasing  the 
sampling  frequency.  For  example,  sampling  at  51 2  Hz  decreases  the 
residue  errors  for  H^z)  to  -0.00211  =fj0.00379. 


(a)  Bode  magnitude  plot  for  H2(s)  (6)  Bode  magnitude  plot  for  H2*(s) 


Fig.  7.  Bode  magnitude  plots  for  Example  2:  (a)  original  system,  and  (b)  the  system  after  removing  the  pole  at  the  origin. 


(a)  Continuous-time  impulse  responses  for  H2*(s)  (b) 


Continuous-time  step  responses  for  H2*(s) 


Fig.  8.  Approximation  of  the  continuous  time  responses  for  Example  2  using  the  DRA. 
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Approximate  discrete-time  impulse  resp.  for  H2*(s) 
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Fig.  9.  System  impulse  response  at  a  sampling  period  of  0.1  s  for  Example  2. 
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Singular  values  of  Hankel  matrix  for  H2*(s) 
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Fig.  10.  Singular  values  of  Hankel  matrix  in  Example  2. 


3.2.  Example  2:  rational  polynomial  transfer  function  with  a  pole 
at  the  origin 

In  this  example,  we  demonstrate  how  to  handle  a  single  pole  at 
the  origin.  The  continuous-time  transfer  function  is  given  by 

H2(s)  =  -(  — l - -)  .  (15) 

s  Vs2  +  6s  +  87 

This  system  has  real  poles  at  0,  2  and  4 rads-1.  We  require  a 
discrete-time  transfer  function  with  a  sampling  period  of  Ts  =  0.1  s. 
Prior  to  Step  1  we  remove  the  pole  at  the  origin.  This  is  accom¬ 
plished  by  first  calculating  the  residue  for  this  pole.  In  this 
polynomial  example,  the  residue  can  be  computed  analytically  as 


Final  model  discrete-time  impulse  response  for  H2(s) 

x  io'3 


2. 


res0  =  lim  sH(s)  =  0.125. 


In  general,  we  find  this  residue  by  selecting  a  very  small  value  for 
s  and  numerically  computing  reso.  The  reduced  transfer  function, 
H*(s)  with  the  pole  at  the  origin  removed  is 


H?(S)  =  - 


If _ I _ I 

s  V  s2  +  6s  +  8  ) 


0.125 

s 


Step  4:  The  state-space  representation  found  in  Step  3  is  aug¬ 
mented  to  include  the  pole  at  the  origin  according  to  Eqs.  (12)  and 
(13).  The  discrete-time  realization  of  H2(s)  is 
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Fig.  7  shows  the  magnitude  plot  of  the  original  system  and  the 
system  with  the  pole  at  the  origin  removed. 

Step  1 :  H*(s)  is  sampled  at  256  Hz  which  is  significantly  greater 
than  the  system  bandwidth.  The  approximate  continuous-time 
impulse  response  is  computed  and  plotted  in  Fig.  8(a).  The  length 
of  the  impulse  response  is  truncated  to  64  samples  which  allows 
for  a  Hankel  matrix  up  to  32  x  32. 

Step  2:  the  approximation  to  the  continuous-time  step  response 
of  H*(s)  is  calculated  as  in  the  first  example  and  plotted  in  Fig.  8(b). 
This  step  response  is  sampled  at  Ts  =  0.1  s,  and  differenced  to  yield 
the  discrete-time  impulse  response,  plotted  in  Fig.  9. 

Step  3 :  The  system  Hankel  matrix  is  generated  from  the  discrete¬ 
time  impulse  response  found  in  Step  2.  Fig.  10  depicts  the  32 
singular  values  of  the  system  Hankel  matrix.  The  first  two  singular 
values  are  two  orders  of  magnitude  greater  than  the  third,  indicat¬ 
ing  that  H*(s)  is  a  second  order  system.  The  Ho-Kalman  algorithm 
generates  the  A ,  B ,  and  C  matrices  after  truncating  all  but  the  first 
two  states.  The  value  of  D  in  this  example  is  0. 


B 


aug  — 


0.1162 

-0.03401 


r  - 

'-'aug 


0.125 


-0.1162  -0.03401 


D  =  [0]. 


Fig.  1 1  shows  close  comparison  of  the  impulse  response  found  from 
the  DRA  and  the  exact  solution.  The  poles  of  the  exact  discrete-time 
system  are  found  at  1 , 0.81 9  and  0.670  with  corresponding  residues 
of  0.125,  -0.0227  and  0.0103.  The  difference  in  residues  between 
the  exact  solution  and  the  DRA  realization  are  0,  -8.88  x  10-5  and 
8.03  xlO-5. 


3.3.  Example  3:  transcendental  transfer  function 

In  the  first  two  examples,  we  used  rational  polynomials  to  illus¬ 
trate  the  DRA  method  where  order  of  the  system  is  known  a  priori , 
and  the  exact  answer  could  be  calculated  analytically.  We  will  now 
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Bode  magnitude  plot  for  H3(s) 


Bode  magnitude  plot  for  H3*(s)  [exact] 


Fig.  12.  Bode  magnitude  responses  Example  3a:  (a)  original  equation;  (b)  with  pole  at  origin  removed. 


(a) 


Continuous-time  impulse  response  for  H3*(s) 


(b) 


Continuous-time  step  response  for  H3*(s) 


Fig.  13.  Continuous  time  responses  for  Example  3a. 


demonstrate  the  DRA  with  an  infinite-order  distributed-parameter 
system.  Specifically  the  algorithm  is  used  on  a  model  of  lithium  dif¬ 
fusion  in  a  spherical  solid  particle  of  an  electrochemical  cell  which 
is  described  by, 

dcs  Ds  3  (  2  9cs  \ 

3 1  r2  dr  y  3r  J  ’ 


with  boundary  conditions, 


-Ds 


3  cs 
dr 
3  cs 
dr 


r= 0 


r=Rs 


=  0 

_  jU(t ) 
asF  ' 


In  these  equations,  cs  is  the  concentration  of  lithium  in  the  elec¬ 
trode  particle  and  r  is  the  radial  distance  from  the  center  of  the 
particle.  The  range  of  r  is  from  0,  the  center  of  the  particle,  to  Rs,  the 
radius  of  solid  electrode  particle.  Ds  is  the  solid  diffusivity  and  F  is 
Faraday’s  constant.  The  specific  interfacial  surface  area  is  as  and/1 
is  the  reaction  current.  The  continuous-time  transfer  function  for 
this  system  was  derived  by  Jacobsen  and  West  [30] 


H3(s)  = 


cs,e(s) 

7W 


Rs  tanh(yft) 
cisDsF  tanh(/3)  -  /3 


(16) 


where  fd  =  Rs  y/s/Ds  and  cs,e(t)  =  cs(Rs,  t)  is  the  surface  concentration 
of  lithium.  The  transfer  function  has  a  pole  at  the  origin.  In  Exam¬ 
ple  3a,  this  pole  is  removed  in  the  same  manner  as  used  by  Smith, 
which  is  to  analytically  subtract  out  the  average  concentration 


[4].  In  Example  3b,  the  pole  is  removed  by  the  numerical  residue 
method.  Ultimately,  both  examples  produce  a  reduced-order, 
discrete-time,  state-space  realization  of  Eq.  (16).  Table  1  lists  the 
parameters  and  values  used  in  this  example. 

3.4.  Example  3a:  diffusion  equation  with  pole  at  origin  removed 
exactly,  analytically 

The  average  concentration  of  lithium  in  the  solid  phase  is  given 
by 


cs,avg(s)  _  ~3  1 

jLi(s)  ~  RsOsFs '  1  J 

Smith  showed  that  subtracting  Eq.  (17)  from  Eq.  (16)  gives  the 
following  transfer  function,  where  Acs,e(s)  =  cs,e(s)  -  cs,aVg(s): 


H|(s)  = 


ACs,e(s) 

jLi(S ) 


Rs  {01  +  3)tanh(/J)  -  3/1 
asDsF  ^2(tanh(/g)  _  p) 


Table  1 

Parameters  used  in  Examples  3a  and  3b. 


Parameter  name 

Interpretation 

Value 

Ts 

Sampling  period 

1  s 

Rs 

Particle  radius 

10-5  m 

Ds 

Diffusivity 

10-12  m2  s_1 

as 

Specific  interfacial  surface  area 

1.74  x  105  m-1 

cs(r ,  0) 

Initial  lithium  concentration 

10,000  mol  ur3 
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Approximate  discrete-time  impulse  resp.  for  H^fs) 


Fig.  14.  Impulse  response  of  the  H*(s )  for  Example  3a. 


Singular  values  of  Hankel  matrix  for  H3*(s) 


Fig.  15.  Singular  values  of  Hankel  matrix  for  Example  3a. 


We  find  a  reduced-order  model  for  H3(s)  by  first  using  the  DRA  on 
H*(s).  The  pole  at  the  origin,  given  by  Eq.  (17),  is  then  added  as  the 
final  step  of  the  procedure,  and  the  final  result  is  a  reduced-order 
state-space  realization  of  H3(s) . 

Step  1 :  The  magnitude  responses  of  H3 (s)  and  H* (s)  are  shown  in 
Fig.  12.  H*(s)  is  sampled  at  256  Hz  for  a  total  of  256  s.  The  approx¬ 
imate  continuous-time  impulse  response  is  shown  in  Fig.  13(a). 
(There  is  no  known  exact  solution  to  the  continuous-time  impulse 
response  for  this  system  against  which  to  compare  this  result.) 

Step  2:  The  approximate  continuous-time  step  response  is  cal¬ 
culated  by  performing  a  cumulative  sum  of  the  impulse  response  of 
Step  1.  Fig.  13(b)  shows  the  approximation  of  the  continuous  time 
impulse  and  step  response.  (Again,  there  is  no  known  exact  solu¬ 
tion  to  the  continuous-time  step  response  for  this  system  against 
which  to  compare  this  result.)  The  approximate  continuous-time 
step  response  is  sampled  at  Ts  =  1  second,  and  differenced  to  pro¬ 
duce  the  discrete  time  impulse  response,  shown  in  Fig.  14. 

Step  3:  The  Hankel  matrix  is  formed,  and  the  singular  values  are 
plotted  in  Fig.  15.  H*(s)  represents  a  distributed-parameter  sys¬ 
tem  that  actually  has  an  infinite  number  of  poles.  However,  we  see 
from  this  plot  that  only  a  few  of  these  poles  are  significant  to  the 
solution.  In  particular,  we  choose  to  use  a  reduced-order  model 
dimension  p  =  2  in  the  results  we  present  here.  This  demonstrates 
a  tradeoff  between  the  complexity  and  accuracy  of  the  solution. 
In  this  example,  the  realization  is  formed  by  truncating  all  but  the 
first  two  states.  Another  approach  would  be  to  reduced  the  order 
to  a  value  higher  than  two  and  do  a  balanced  residualization  [27] 
in  Step  4. 

Step  4:  The  state-space  realization  derived  in  Step  3  is  aug¬ 
mented  with  the  integrator  state  to  give  the  final  third-order 


state-space  model  of  the  diffusion  equation.  The  final  realization 
is  given  by, 


A 


aug  - 


B 


aug  — 


1 

0  0 

0 

0.4695  0.3296 

0 

0.3296  0.4355 

1 


7.094  x  10“3 
-1.698  x  10"3 


6'aug  — 


-1.787  x  10“5 


-7.094  xl0“3  1.698  xl0“3 


,  D  =  [0], 


The  output  of  this  discrete-time  realization  to  a  10  s  discharge  fol¬ 
lowed  by  a  10  s  rest  is  shown  in  Fig.  16(a).  This  is  validated  against 
a  solution  computed  by  a  1-D  parabolic-elliptic  partial  differential 
equation  solver.  Errors  between  the  PDE  solution  and  DRA  solu¬ 
tion  are  shown  in  Fig.  16(b).  The  DRA  third-order  model  accurately 
models  the  system  behavior. 


3.5.  Example  3b:  diffusion  equation  with  pole  at  origin  removed 
numerically 

In  Example  3a,  the  pole  at  the  origin  is  removed  exactly  by  sub¬ 
tracting  out  the  average  concentration.  In  general,  it  may  be  difficult 
or  impossible  to  find  a  closed-form  algebraic  function  that  allows 
the  pole  at  the  origin  to  be  removed  exactly.  For  those  situations, 


Fig.  16.  Simulation  results  of  a  10  s  discharge  followed  by  a  10  s  rest  for  Example  3a. 
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Fig.  17.  Comparison  of  final  results  by  removing  the  pole  at  the  origin  by  analytic  method  and  by  the  residue  method  for  Example  3b. 


the  pole  is  numerically  removed  prior  to  Step  1.  The  residue  for  the 
pole  at  the  origin  is  found  with 


res0  =  lims 

s^O 


Cs,e(s) 

iLKs) 


-3  Rs 
asDsF 


which  we  find  using  Mathematica.  We  also  achieve  similar  results 
using  a  very  small  value  for  s  (e.g.,  s  =  10-10).  For  this  example 
res0  =  - 1.787  x  10-5.  This  result  is  subtracted  from  the  original 
transfer  function  to  give 


H*(s)  = 


Acs,e(s) 
jLi(s ) 

Rs  tanh(/3) 
asDsF  tanh (0)  -  p 


-1.787  x  10"5 
s 


The  same  sampling  frequency  (256  Hz)  and  the  same  impulse- 
response  length  (256  s)  are  used.  The  singular  values  are  very 
similar  to  Example  3a  and  again  a  third-order  state  space  model 
is  generated.  Simulation  results  in  Fig.  17(a)  show  excellent  agree¬ 
ment  between  the  analytic  method  and  the  residue  method  for 
accounting  for  a  pole  at  the  origin.  Figure  17(b)  shows  the  differ¬ 
ence  between  the  Example  3a  and  Example  3b  solutions,  which  are 
on  the  order  of  1 0-3  mol  rrr3. 


4.  Conclusions 

In  this  paper,  we  introduced  the  discrete-time  realization 
algorithm  (DRA)  which  is  used  to  arrive  at  a  discrete-time, 
reduced-order  state-space  model  of  a  Laplace  domain  transfer 
function.  This  method  is  based  on  deriving  an  approximation  to 
the  impulse  response  of  the  system  and  then  using  this  impulse 
response  with  the  Ho-Kalman  algorithm.  Although  the  examples 
were  for  single-input-single-output  systems  this,  method  works 
with  multiple-input  multiple-output  systems  as  well.  Note  that 
the  steps  of  the  DRA  require  only  standard  linear-algebra  and 
signal-processing  routines.  The  method  finds  the  globally  opti¬ 
mal  reduced-order  discrete-time  approximation  to  the  original 
continuous-time  system.  It  does  not  require  nonlinear  optimiza¬ 
tion,  and  does  not  require  iteration.  In  a  future  work  we  will  use 
the  DRA  to  derive  a  reduced-order  single-input  multiple-output 
lithium-ion  cell  model. 
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